#include <iostream>
#include <fstream>
#include <string>
#include <vector>
#include <algorithm>
#include <ctime>
#include <cstdlib>
#include <cstdio>
#include <cmath>
#include <unistd.h>
#include <math.h>
#include <exception>

#include <sstream>
#include <time.h>
#include <sys/time.h>

#include "PNN1Inference.h"

using namespace std;

#define SM (CLS / sizeof (double))

int confusion[30][30];
int confusion_sum[30][30];
int contactStat[30];
int contactStat_sum[30];
int sum[30];

SEQUENCE::SEQUENCE(int length_seq, PNN1Inference* pModel, string name)
{
  m_pModel = pModel;
  p_name= name;
  this->length_seq = length_seq;
  length_contact=length_seq*length_seq;
  dim_features=m_pModel->dim_features;
  window_size= m_pModel->window_size;
  dim_one_pos=m_pModel->dim_one_pos;
  num_gates = m_pModel->num_gates;
  num_states = m_pModel->num_states;
  num_labels=m_pModel->num_labels;
  dim_pairwise_features = m_pModel->dim_pairwise_features;
  gates=NULL;
  arrVi=NULL;
  Partition=NULL;
  obs_feature = new Score*[length_seq];
  for(int i=0;i<length_seq;i++)
    obs_feature[i] = new Score[dim_one_pos];
  s_start=0;
  s_end=length_seq-1;
  output_center=NULL;
  predicted_label=NULL;
  Rg=2.2*pow(length_seq, 0.38);
  pFeatures=new Score[dim_features];
  obs_exist_flag=true;
  m_sumx=0, m_sumx2=0, m_sumy=0, m_sumxy=0, m_sumy2=0, m_sumR=0, m_len=0;
}

SEQUENCE::~SEQUENCE(){
  delete arrVi;
  delete gates;
  delete Partition;
  delete obs_label;

  for(int i=0;i<length_seq;i++)
    delete obs_feature[i];
  delete obs_feature;
  delete output_center;
  delete predicted_label;
  delete pFeatures;
}

Score SEQUENCE::MAP(){ // Posterior Decoding (Marginal Probability Decoder)
  Score bestScore, obj = 0;
  Score *vi2;
  int i, t, bestState=0; 

  CalcPartition();
  for(t=nDataPoints-1; t>=0; t--) {
    //cerr << " MAP: " << t << ": " << obs_label[t] << endl;
    obj+= (*arrVi)(t, obs_label[t]) - Partition[t];
  }

  for(t=nDataPoints-1, vi2=&((*arrVi)(t, 0)); t>=0; t--, vi2-=arrVi->cols){
    bestState=0; 
    bestScore=vi2[0];
    for(i=1; i< num_states; i++){
      if (vi2[i]>bestScore){
        bestScore=vi2[i];
        bestState=i;
      }
    }
    predicted_label[t]=bestState;
  }
  return obj;
}

void SEQUENCE::CalcPartition()
{
  int t, i, t2, i2, lastt=nDataPoints/SM*SM, lasti=num_states/SM*SM;
  Score score, *vi2, *P2;

  if (Partition==NULL) 
    Partition = new Score[nDataPoints];
  for(t=0;t<nDataPoints;t++)
    Partition[t] = (Score)LogScore_ZERO;
// LogScore_PLUS_EQUALS(Partition[t], (*arrVi)(t, state));
  for(t=0;t<nDataPoints;t+=SM)
    for(i=0; i< num_states; i+=SM)
      for ( t2 = 0, vi2=&((*arrVi)(t, i)), P2=&Partition[t]; 
            t2 < (t<lastt ? SM : nDataPoints%SM);
            t2++, vi2 += arrVi->cols)
        for ( i2 = 0; i2 < (i<lasti ? SM : num_states%SM); i2++){
          score = vi2[i2];
          LogScore_PLUS_EQUALS(P2[t2], score);
          //cerr << "  Partition(" << t << "," << state << ")=" << score << " : " << Partition[t] << endl;
        }
  //cerr << endl;
}

void SEQUENCE::ComputeVi()
{
  int pos, pos2, l, l2, g, g2,
      lastpos=nDataPoints/SM*SM, lastl=num_states/SM*SM, lastg=num_gates/SM*SM;
  Score *gates2, *vi2, *pW2;

  if (arrVi==NULL)
    arrVi = new ScoreMatrix(nDataPoints, num_states, "arrVi");
  arrVi->resize(nDataPoints, num_states); //arrVi->Fill(LogScore_ZERO);

// (*arrVi)(pos,l) += m_pModel->weights[GetWeightLAddr(l, g)]*(*gates)(pos,g);

  for ( pos = 0; pos<nDataPoints; pos+=SM)
  for ( l = 0; l<num_states; l+=SM)
  for ( g = 0; g<num_gates; g+=SM) {
    for ( pos2 = 0, gates2 = &((*gates)(pos, g)), vi2 = &((*arrVi)(pos, l));
          pos2 < (pos<lastpos ? SM : nDataPoints%SM); 
          pos2++, vi2 += arrVi->cols, gates2 += gates->cols)
      for ( l2 = 0, pW2 = m_pModel->GetWeightLAddr(l,g); 
            l2 < (l<lastl ? SM : num_states%SM); 
            l2++, pW2+=num_gates)
      {
        for ( g2 = 0; g2 < (g<lastg ? SM : num_gates%SM); g2++)
        {
          vi2[l2] += pW2[g2] * gates2[g2];
          // if (vi2[l2]>10000)
          //   cerr << " ComputeVi: Vi("<<pos+pos2<<", "<<l+l2<<")= W("<<l+l2<<","
          //        << g+g2 << ") * gates("<<pos+pos2<<","<<g+g2<<"):: "<<vi2[l2]
          //        << " += " << pW2[g2] << " * "<<gates2[g2]<<endl; 
        } // exit(0);
      }
  }
  
  //for (pos=0; pos<length_contact; pos++)
  //  for (l=0; l<num_states; l++)
  //    cerr << "  ComputeVi(" << pos << ", " << l << ")=" << (*arrVi)( pos, l) << endl;
}

void SEQUENCE::Predict() 
{
  ComputeGates();
  ComputeVi();
  //cerr << "  Predictions " << nDataPoints << endl;
  MAP();

  for (int t=0; t<nDataPoints; t++) {
    cout << predicted_label[t] << " " ;
    for (int l=0; l<num_labels; l++)
      cout << exp((*arrVi)(t, l) - Partition[t]) << " ";
    cout << endl;
  }
}

SingletonSequence::SingletonSequence(int len, PNN1Inference* pModel, string name)
  : SEQUENCE(len, pModel, name)
{
  obs_label = new int[length_seq];
  _features=NULL;
}

SingletonSequence::~SingletonSequence()
{
  delete _features;
}

int SingletonSequence::makeFeatures(){
  //build features for local windows  
  _features = new ScoreMatrix(length_seq, dim_features, "_features");
  for(int t=0;t<length_seq;t++){
    if (t<window_size/2)
      (*_features)(t, t) = 1;
    else if (length_seq-1-t<window_size/2)
      (*_features)(t, length_seq-1-t) = 1;

    int pivot = window_size/2*2;
    for(int i=0; i< window_size; i++){
      int offset = t + i - window_size/2;
      if(offset <0 or offset >=length_seq){
        for(int j=0; j<dim_one_pos; j++)
          (*_features)(t, pivot) = 0, pivot++;
      } else {
        for(int j=0; j<dim_one_pos; j++)
          (*_features)(t, pivot) = obs_feature[offset][j], pivot++;
      }
    }
    (*_features)(t, pivot) = 1; // last one for the bias term
    // cerr << " ### last pivot=" << pivot << endl;
  }

  nDataPoints = length_seq;
  predicted_label = new int[nDataPoints];
  return 0;
}

void SingletonSequence::ComputeGates()
{
  Score *pfb2, *ocenter, *pW2;
  int t, t2, g, g2, k, k2,
      lasttseq=length_seq/SM*SM, lastg=num_gates/SM*SM, lastk=dim_features/SM*SM;

  struct timeval ts, te;
  gettimeofday(&ts, 0);

  if (gates==NULL){
    gates = new ScoreMatrix(length_seq, num_gates, "gates");
    output_center = new ScoreMatrix(num_gates,length_seq, "output_center");
  } else {
    output_center->resize(num_gates, length_seq); //output_center->Fill(0);
  }
  for (t=0; t<length_seq; t+=SM)
  for (g=0; g<num_gates; g+=SM)
  for (k=0; k<dim_features; k+=SM) 
    for ( g2 = 0, pW2 = m_pModel->GetWeightGAddr(g,k), ocenter = &((*output_center)(g, t));
          g2 < (g<lastg ? SM : num_gates%SM);
          g2++, pW2 += dim_features, ocenter += output_center->cols)
    {
      for ( t2 = 0, pfb2 = &((*_features)(t, k)); 
            t2 < (t<lasttseq ? SM : length_seq%SM);
            t2++, pfb2 += _features->cols) 
      {
        for ( k2 = 0; k2 < (k<lastk ? SM : dim_features%SM); k2++) 
        {
          ocenter[t2] += pW2[k2] * pfb2[k2];
        }
      } //cerr << "  GetGate(" << gate << ", " << pos << ")=" << output << endl; //exit(0);
    }
  for ( t=0; t<length_seq; t++)
    for ( g=0; g<num_gates; g++){
      (*gates)(t, g) = Gate((*output_center)(g, t));
      // if (t==30) cerr << "   gate("<< t <<"," << g <<")= " << (*gates)(t, g) << endl;
    }
  gettimeofday(&te, 0);
  double dif = difftime (te.tv_sec,ts.tv_sec);
  //cerr << "Gates: " << dif*1000+te.tv_usec/1000-ts.tv_usec/1000 << " miliseconds." << endl;
}

void SingletonSequence::ComputeTestAccuracy(){
  m_pModel->totalPos += nDataPoints;
  // comparison
  for(int t=0; t<nDataPoints;t++) {
    if(obs_label[t]==predicted_label[t])
      m_pModel->totalCorrect++;
    confusion[predicted_label[t]][obs_label[t]]++;
  }
  // for ( int t=0; t<nDataPoints;t++)
  //   cerr << predicted_label[t] << ":" << obs_label[t] << " "; 
  // cerr << "=="<<m_pModel->totalCorrect << endl;
}

PairwiseSequence::PairwiseSequence(int len, PNN1Inference* pModel, string name)
  : SEQUENCE(len, pModel, name)
{
  contact_DIST=7;
  IdxMap=new int[length_seq*length_seq];
  memset(IdxMap, -1, sizeof(int)*length_seq*length_seq);
  rIdxMap=new int[length_seq*length_seq];
  centerMap=new int[length_seq*length_seq];
  IdxMapCenter=new int[length_seq];
  IdxMapCofactor=new int[length_seq];
  for (int i=0; i<length_seq; i++)
    IdxMapCofactor[i]=IdxMapCenter[i]=-1;
  output_cofactor=NULL;
  distances=new ScoreMatrix(length_seq, length_seq, "distances");
  distances->Fill(-1);
  dist_pred=new ScoreMatrix(length_seq, length_seq, "dist_pred");
  obs_label=NULL;
  pairwise_features=NULL;
  obs_pairwise_features=NULL;
  if (dim_pairwise_features)
    obs_pairwise_features=new Array3D<double>(length_seq, length_seq, dim_pairwise_features);
}

PairwiseSequence::~PairwiseSequence()
{
  delete output_cofactor;

  delete IdxMap; 
  delete rIdxMap;
  delete centerMap;
  delete IdxMapCenter;
  delete IdxMapCofactor;
  delete distances;
  delete obs_pairwise_features;
  delete pairwise_features;
  delete dist_pred;
}

int PairwiseSequence::makeFeatures(){
  if (obs_label==NULL) 
    obs_label = new int[length_contact];
  int df = length_contact*2;
  //cerr << proc_id << "@ df=" << df << ", len_contact=" << length_contact << ", gaps_size=" << gaps.size();
  _features = new int[df];
  
  int totalContact=0;

  //cerr << " window_size=" << window_size << ", len_seq=" << length_seq << endl;
  //build features for local windows  
  for(int t=0;t<length_seq-contact_DIST;t++){
      bool bGap=false;
      for (int g=0; g<(int)gaps.size(); g++){
        if (abs(gaps[g]-t)<=window_size/2){
      bGap=true; break;
    }
      }
      if(bGap)
    continue;
      //for(int u=t+7;u<length_seq-window_size/2;u++){
      for(int u=t+contact_DIST;u<length_seq;u++){
    if (totalContact>=length_contact-1)
      break;
    bGap=false;
    for (int g=0; g<(int)gaps.size(); g++){
      if (abs(gaps[g]-u)<=window_size/2){
        bGap=true; break;
      }
    }
    if(bGap) 
      continue;
    int label=(int)(*distances)(t,u);
    if (label<0)
      continue;
    obs_label[totalContact]= label;
    // if (!proc_id)
    //   cerr << "    Adding contact " << totalContact << ":" << p_name << " (" << t << ", " << u 
    //        << ") with dist=" << (*distances)(t,u)<< " labeled " << obs_label[totalContact] << endl;
    IdxMap[u*length_seq+t]=IdxMap[t*length_seq+u]=totalContact;
    centerMap[totalContact]=t;
    IdxMapCenter[t]=totalContact;
    rIdxMap[totalContact]=u;
    IdxMapCofactor[u]=totalContact;
    int pivot = totalContact*2;
    totalContact++;
    int offset1 = t-window_size/2;
    int offset2 = u-window_size/2;
    _features[pivot++]=offset1;
    _features[pivot++]=offset2;
      }
  }
  delete distances; distances=NULL;
  nDataPoints=length_contact=totalContact; 
  if (length_contact==0)
     return -1; //cerr << proc_id << " Making features..." << p_name<< endl;
  predicted_label = new int[nDataPoints];
  pairwise_features = new ScoreMatrix(dim_pairwise_features+1, nDataPoints, "pairwise_features");
  for(int t=0;t<length_contact; t++){
    int offset1=centerMap[t], offset2=rIdxMap[t];
    (*pairwise_features)(0,t)= (offset2-offset1+0.000001)/length_seq;
    for (int d=1; d<=dim_pairwise_features; d++)
      (*pairwise_features)(d,t)=(*obs_pairwise_features)(offset1, offset2, d-1);
  }
  delete obs_pairwise_features;
  obs_pairwise_features = NULL;
  return 0;
}

// Features: ends vector(2*[window_size/2]), features read from input 
// file (dim_one_pos*window_size), sequence distance; sequence distance/length_seq.
Score PairwiseSequence::getFeatures(int pos, int i){
  int pivot=pos*2;
  int offset1 = _features[pivot];
  int offset2 = _features[pivot+1];
  int iif=i-window_size/2*2;
  int idx = iif / m_pModel->dim_one_pos;
  Score ret=1; // when i==dim_features-3-dim_pairwise_features
  if (i<window_size/2*2){
    offset1+=window_size/2;
    offset2+=window_size/2;
    //if (proc_id==0) cerr << "offset1="<<offset1<<", offset2="<<offset2<<endl;
    if (offset1<window_size/2 and i==offset1)
      return 1;
    else if (length_seq-offset2==window_size/2*2-i and offset2>=length_seq-window_size/2){
      //if (proc_id==0) cerr << length_seq << "-" << offset2 << "="<<window_size/2*2<<"-" <<i << endl; 
      return 1;
    } else 
      return 0;
  } else if (idx<window_size){
    if (offset1+idx<0 or offset1+idx>=length_seq) return 0; // training Secondary Structure: allow termial residues
    ret= obs_feature[offset1+idx][iif%m_pModel->dim_one_pos];
  } else if (i<dim_features-3-dim_pairwise_features){
    idx-=m_pModel->window_size;
    //cerr << idx << " - " << i << "*"<<dim_features <<" : " << i%m_pModel->dim_one_pos<<" : " << length_seq << " : " << offset1 << "  : " << offset2;
    if (offset2+idx<0 or offset2+idx>=length_seq) return 0; // training Secondary Structure: allow termial residues
    ret= obs_feature[offset2+idx][iif%m_pModel->dim_one_pos];
    //cerr << " : " << ret << endl;
  } else if (i==dim_features-2-dim_pairwise_features)
    ret = Rg/10;
  else if (i==dim_features-1-dim_pairwise_features)
    ret = (offset2-offset1+0.000001)/length_seq;

  return ret;
}

void PairwiseSequence::ComputeGates()
{
  ScoreMatrix* feature_buffer = new ScoreMatrix(dim_features, length_seq, "feature_buffer");
  ScoreMatrix* feature_buffer_cofactor = new ScoreMatrix(dim_features, length_seq, "feature_buffer_cofactor");
  Score *pfb2, *pfbc2, *ocenter, *ocofactor;
  //int pos;
  int gate, p, t, t2, g, g2, k, k2, kk, 
      lastt = length_contact/SM*SM, lasttseq=length_seq/SM*SM, 
      lastg = num_gates/SM*SM, lastk=(dim_features-1-dim_pairwise_features)/SM*SM;
  Score* gates2, *pW2, output;

  struct timeval ts, te;
  gettimeofday(&ts, 0);

  if (gates==NULL){
    gates = new ScoreMatrix(length_contact, num_gates, "gates");
    output_center = new ScoreMatrix(num_gates,length_seq, "output_center");
    output_cofactor = new ScoreMatrix(num_gates,length_seq, "output_cofactor");
  } else {
    // gates->resize(length_contact, num_gates); //gates->Fill(LogScore_ZERO);
    output_center->resize(num_gates, length_seq); //output_center->Fill(0);
    output_cofactor->resize(num_gates, length_seq); // output_cofactor->Fill(0);
  }
  for(t=0; t<length_seq; t++){
    int idx=IdxMapCenter[t];
    int coidx=IdxMapCofactor[t];
    for(k=0;k<dim_features-1;k++){
      if (k<window_size/2 or (k>=window_size/2*2 and k<dim_features/2+window_size/2-2)){
        if (idx<0) { continue;}
        (*feature_buffer)(k, t)= getFeatures(idx, k);
      } else {
        if (coidx<0) continue;
        (*feature_buffer_cofactor)(k, t) = getFeatures(coidx, k);
      }
    }
  }
  for ( t = 0; t<length_seq; t+=SM) 
  for ( g = 0; g<num_gates; g+=SM)
  for ( k = 0; k<dim_features-1; k+=SM) 
    for ( g2 = 0, pW2 = m_pModel->GetWeightGAddr(g,k),
                  ocenter=&((*output_center)(g, t)), 
                  ocofactor=&((*output_cofactor)(g, t));
          g2 < (g<lastg ? SM : num_gates%SM); g2++)
    {
      for ( k2 = 0, pfb2 = &((*feature_buffer)(k, t)),
                    pfbc2=&((*feature_buffer_cofactor)(k, t));
            k2 < (k<lastk ? SM : (dim_features-1-dim_pairwise_features)%SM);
            k2++ ) 
      {
        kk = k + k2;
        if (kk<window_size/2 or (kk>=window_size/2*2 and kk<dim_features/2+window_size/2-2)){
          for(t2=0; t2<(t<lasttseq?SM:length_seq%SM); t2++) 
            ocenter[t2] += pW2[k2] * pfb2[t2];
        } else {
          for(t2=0; t2<(t<lasttseq?SM:length_seq%SM); t2++)
            ocofactor[t2] +=pW2[k2] * pfbc2[t2];
        }
        pfb2+=feature_buffer->cols;
        pfbc2+=feature_buffer_cofactor->cols;
      } //cerr << "  GetGate(" << gate << ", " << pos << ")=" << output << endl; //exit(0);
      pW2+=dim_features;
      ocenter+=output_center->cols;
      ocofactor+=output_cofactor->cols;
    }
  for (t=0; t<length_contact; t+=SM)
  for (g=0; g<num_gates; g+=SM)
    for (t2=0, gates2=&((*gates)(t, g)); t2<(t<lastt?SM:length_contact%SM); t2++, gates2+=gates->cols){
      int center=centerMap[t+t2];
      int cofactor=rIdxMap[t+t2];
       for (g2=0; g2<(g<lastg?SM:num_gates%SM); g2++) {
      gate=g+g2; //pos=t+t2;
      pW2=m_pModel->GetWeightGAddr(gate,0);

      output=(*output_center)(gate, center) + (*output_cofactor)(gate, cofactor);
      for (p=0; p<dim_pairwise_features+1; p++)
        output+=pW2[dim_features-dim_pairwise_features-1+p]*(*pairwise_features)(p,t);

      gates2[g2] = Gate(output);
      //if (pos=30) cerr << "   gate("<< gate <<"," << pos <<")= " << gates2[g2] << " " << output << endl;
       }
    }
  gettimeofday(&te, 0);
  double dif = difftime (te.tv_sec,ts.tv_sec);
  //cerr << "Gates: " << dif*1000+te.tv_usec/1000-ts.tv_usec/1000 << " miliseconds." << endl;
  delete feature_buffer;
  delete feature_buffer_cofactor;
}

void PairwiseSequence::ComputeTestAccuracy(){
  m_pModel->totalPos += length_contact;
  for(int t=0; t<length_contact;t++) {
    if(obs_label[t]==predicted_label[t])
      m_pModel->totalCorrect++;
    confusion[predicted_label[t]][obs_label[t]]++;
  }
}

/*
// func=0: contact prediction for top L/5, 1-top L/10, 2-top 5; 3-energy function for FM, 4-energy function for TB; 5-contact number
void PairwiseSequence::Predict() 
{
  ComputeGates();
  ComputeVi();
  CalcPartition();
  int pred_positive=0, top=0, TOP_most=2;
  if (m_pModel->startRes<0)
    m_pModel->startRes=1;
  if (m_pModel->endRes<1)
    m_pModel->endRes=length_seq;
  int totalRes=min(m_pModel->endRes, length_seq) -m_pModel->startRes+1;
  int needed=length_contact, nonContact=-1;;
  while (pred_positive<needed and top<TOP_most)
    pred_positive=MAP(top++, needed, nonContact);

  for(int i=0;i<num_labels;i++) for(int j=0;j<num_labels;j++) { confusion[i][j]=0; confusionEnergy[i][j]=0;}

  // comparison
  for(int t=0; t<length_contact;t++) {
    int x=centerMap[t], y=rIdxMap[t];
    int label=predicted_label[t];
    (*dist_pred)(x,y)=(*dist_pred)(y,x)=label;
    //cerr << t << "(" << x <<","<<y<< "): " << predicted_label[t] << endl;
  }
}
*/
/////////////////////  implementation of the class ContactNumberSequence ///////////////////////////////////////

DistancePotentialSequence::DistancePotentialSequence(int len, PNN1Inference* pModel, string name) 
  : PairwiseSequence(len, pModel, name)
{
}

int DistancePotentialSequence::makeFeatures(){
  if (obs_label==NULL) 
    obs_label = new int[length_contact];
  int df = length_contact*2;
  //cerr << proc_id << "@ df=" << df << ", len_contact=" << length_contact << ", gaps_size=" << gaps.size();
  _features = new int[df];
  
  int totalContact=0;

  //cerr << " window_size=" << window_size << ", len_seq=" << length_seq << endl;
  //build features for local windows  
  int negativeLabel=0;
  //for(int t=max(s_start, window_size/2);t<min(s_end, length_seq-window_size/2);t++){
  //for(int t=window_size/2;t<length_seq-window_size/2;t++){
  for(int t=0;t<length_seq-contact_DIST;t++){
      bool bGap=false;
      for (int g=0; g<(int)gaps.size(); g++){
        if (abs(gaps[g]-t)<=window_size/2){
      bGap=true; break;
    }
      }
      if(bGap)
    continue;
      //for(int u=t+7;u<length_seq-window_size/2;u++){
      for(int u=t+contact_DIST;u<length_seq;u++){
    if (totalContact>=length_contact-1)
      break;
    bGap=false;
    for (int g=0; g<(int)gaps.size(); g++){
      if (abs(gaps[g]-u)<=window_size/2){
        bGap=true; break;
      }
    }
    if(bGap) 
      continue;
    int label = (int)(((*distances)(t,u)-3)*1);
    //int label1 = (int)(((*distances)(t,u)-7)*1), max_label=12;   // 8 states: 3-7.99, 8-8.99, ..., 13-13.99, 14-inf
    //int label1 = (int)(((*distances)(t,u)-6)*1), max_label=12; // 9 states: 3-6.99, 7-8.99, 8-8.99, ..., 13-13.99, 14-inf
    //int label1 = (int)(((*distances)(t,u)-5)*1), max_label=12; // 10 states: 3-5.99, 6-6.99, 7-7.99, 8-8.99, ..., 13-13.99, 14-inf
    //int label1 = (int)(((*distances)(t,u)-3)*1), max_label=12; // 12 states: 3-3.99, 4-4.99,5-5.99, 6-6.99, 7-7.99, 8-8.99, ..., 13-13.99, 14-inf
    int label1 = (int)(((*distances)(t,u)-4)*1), max_label=13; // 12 states: 3-4.99, 5-5.99, 6-6.99, 7-7.99, 8-8.99, ..., 13-13.99, 14-14.99, 15-inf
    //int label1 = (int)(((*distances)(t,u)-4)*1), max_label=18; // 17 states: 3-4.99, 5-5.99, 6-6.99, ..., 18-18.99, 19-19.99, 20-inf
    if (label1>0 and label<max_label)
      obs_label[totalContact]=label1;
      //obs_label[totalContact]=(label1+1)/2;
    else if (label>=0 and label<max_label)
      obs_label[totalContact]=0;
    else { // when label is negative, it means distance is over 21A for 17 states or over 16A for 12 states, so it's obs_label is the last state possible
      obs_label[totalContact]=num_states-1;
      negativeLabel++;
    } // if (!proc_id) cerr << "    Adding contact " << totalContact << ":" << p_name << " (" << t << ", " << u << ") with dist=" << (*distances)(t,u)<< " labeled " << obs_label[totalContact] << endl;
    IdxMap[u*length_seq+t]=IdxMap[t*length_seq+u]=totalContact;
    centerMap[totalContact]=t;
    IdxMapCenter[t]=totalContact;
    rIdxMap[totalContact]=u;
    IdxMapCofactor[u]=totalContact;
    int pivot = totalContact*2;
    totalContact++;
    int offset1 = t-window_size/2;
    int offset2 = u-window_size/2;
    _features[pivot++]=offset1;
    _features[pivot++]=offset2;
      }
  }
  //cerr << proc_id << "@  Ratio: " << negativeLabel << "/" << totalContact << endl;
  //if (length_contact <=totalContact+2)
  //  cerr << "  ^^^ Got more contact than expected: " << length_contact << " : " << totalContact<< endl;
  delete distances; distances=NULL;
  nDataPoints=length_contact=totalContact; 
  if (length_contact==0)
     return -1; //cerr << proc_id << " Making features..." << p_name<< endl;
  predicted_label = new int[nDataPoints];
  pairwise_features = new ScoreMatrix(1, nDataPoints, "pairwise_features");
  for(int t=0;t<length_contact; t++){
    int offset1=centerMap[t], offset2=rIdxMap[t];
    (*pairwise_features)(0,t)= (offset2-offset1+0.000001)/length_seq;
  }

  return 0;
}

Score DistancePotentialSequence::getFeatures(int pos, int i){
  int pivot=pos*2;
  int offset1 = _features[pivot];
  int offset2 = _features[pivot+1];
  int iif=i-window_size/2*2;
  int idx = iif / m_pModel->dim_one_pos;
  Score ret=1; // when i==dim_features-3
  if (i<window_size/2*2){
    offset1+=window_size/2;
    offset2+=window_size/2;
    //if (proc_id==0) cerr << "offset1="<<offset1<<", offset2="<<offset2<<endl;
    if (offset1<window_size/2 and i==offset1)
      return 1;
    else if (length_seq-offset2==window_size/2*2-i and offset2>=length_seq-window_size/2){
      //if (proc_id==0) cerr << length_seq << "-" << offset2 << "="<<window_size/2*2<<"-" <<i << endl; 
      return 1;
    } else 
      return 0;
  } else if (idx<window_size){
    if (offset1+idx<0 or offset1+idx>=length_seq) return 0; // training Secondary Structure: allow termial residues
    ret= obs_feature[offset1+idx][iif%m_pModel->dim_one_pos];
  } else if (i<dim_features-3){
    idx-=m_pModel->window_size;
    //cerr << idx << " - " << i << "*"<<dim_features <<" : " << i%m_pModel->dim_one_pos<<" : " << length_seq << " : " << offset1 << "  : " << offset2;
    if (offset2+idx<0 or offset2+idx>=length_seq) return 0; // training Secondary Structure: allow termial residues
    ret= obs_feature[offset2+idx][iif%m_pModel->dim_one_pos];
    //cerr << " : " << ret << endl;
  } else if (i==dim_features-2)
    ret = Rg/10;
  else if (i==dim_features-1)
    ret = (offset2-offset1+0.000001)/length_seq;

  return ret;
}

/////////////////////  implementation of the class ContactNumberSequence ///////////////////////////////////////

ContactNumberSequence::ContactNumberSequence(int len, PNN1Inference* pModel, string name)
  : PairwiseSequence(len, pModel, name)
{
  contact_DIST=3;
  obs_cnum=new int[len];
  memset(obs_cnum, 0, sizeof(int)*len);
  pred_cnum = new int[length_seq];
  memset(pred_cnum, 0, sizeof(int)*length_seq);
  logReferenceState = new double[_MAX_CONTACT];
  double refProbs[]={2.66, 3.33, 9.41, 10.41, 14.17, 12.76, 11.23, 10.09, 9.33, 16.61};
  for (int i=0; i<_MAX_CONTACT; i++)
    logReferenceState[i]=log(refProbs[i]/100);
  prob_residue2state=new ScoreMatrix(length_seq, _MAX_CONTACT, "residue2state");
  log_prob=new ScoreMatrix(length_seq, _MAX_CONTACT, "log_prob");
}

ContactNumberSequence::~ContactNumberSequence()
{
  delete obs_cnum;
  delete pred_cnum;
  delete prob_residue2state;
  delete log_prob;
}

int ContactNumberSequence::makeFeatures(){
  if (obs_label==NULL) 
    obs_label = new int[length_contact];
  int df = length_contact*2;
  // cerr << proc_id << "@ df=" << df << ", len_contact=" << length_contact << ", gaps_size=" << gaps.size();
  _features = new int[df];
  
  int totalContact=0;
  int window_size = m_pModel->window_size;

  //cerr << " window_size=" << window_size << ", len_seq=" << length_seq << endl;
  int* nContacts=new int[length_seq];
  int* nContactsR=new int[length_seq];
  for(int t=0;t<length_seq-contact_DIST;t++){
    bool bGap=false;
    for (int g=0; g<(int)gaps.size(); g++){
      if (abs(gaps[g]-t)<=window_size/2){
        bGap=true; break;
      }
    }
    if(bGap)
      continue;

    memset(nContacts, 0, sizeof(int)*length_seq);
    memset(nContactsR, 0, sizeof(int)*length_seq);

    int nn=0;
    for ( int u=t+contact_DIST;u<length_seq;u++){
      nContacts[u]=nContacts[u-1];
      bGap=false;
      for (int g=0; g<(int)gaps.size(); g++)
        if (abs(gaps[g]-u)<=window_size/2) { 
          bGap=true; break; 
        }

      if (bGap) continue;

// For predicting the contact number between residue t and all the residues from t+contact_DIST to u
      if ((*distances)(t,u)<=8.0 and (*distances)(t,u)>2)
        nContacts[u]++, nn++;
      if (nContacts[u]<_MAX_CONTACT-1)
        obs_label[totalContact] = _MAX_CONTACT-1-nContacts[u];
      else obs_label[totalContact] = 0; // merge those bins with contact number more than 9 as one bin
      //if (proc_id==0) cerr << "(" << t << ", " << u << "): " << (*distances)(t,u) << "  " << nContacts[u] << endl;
      contactStat[obs_label[totalContact]]++;

      IdxMap[t*length_seq+u] = totalContact;
      centerMap[totalContact] = t; IdxMapCenter[t]=totalContact;
      rIdxMap[totalContact]=u; IdxMapCofactor[u]=totalContact;

      int pivot = totalContact*2;
      totalContact++;

      int offset1 = t-m_pModel->window_size/2, offset2 = u-m_pModel->window_size/2;
      _features[pivot++]=offset1; _features[pivot++]=offset2;
    }
    for(int u=t-contact_DIST;u>=0;u--){
      nContactsR[u]=nContactsR[u+1];
      bGap=false;
      for (int g=0; g<(int)gaps.size(); g++)
        if (abs(gaps[g]-u)<=window_size/2){ bGap=true; break; }
      if(bGap) continue;
// For predicting the contact number between residue t and all the residues from t+contact_DIST to u
      if ((*distances)(t,u)<=8.0 and (*distances)(t,u)>2)
        nContactsR[u]++, nn++;
      if (nContactsR[u]<_MAX_CONTACT-1)
        obs_label[totalContact] = _MAX_CONTACT-1-nContactsR[u];
      else
        obs_label[totalContact]=0; // merge those bins with contact number more than 9 as one bin
      //if (proc_id==0) cerr << "(" << t << ", " << u << "): " << (*distances)(t,u) << "  " << nContacts[u] << endl;
      contactStat[obs_label[totalContact]]++;

      IdxMap[t*length_seq+u]=totalContact;
      centerMap[totalContact]=t; IdxMapCenter[t]=totalContact; 
      rIdxMap[totalContact]=u; IdxMapCofactor[u]=totalContact;

      int pivot = totalContact*2;
      totalContact++;
      int offset1 = t-m_pModel->window_size/2, offset2 = u-m_pModel->window_size/2;
      _features[pivot++]=offset1; _features[pivot++]=offset2;
    }
    nn=nContactsR[0]+nContacts[length_seq-1];
    obs_cnum[t]=min(_MAX_CONTACT-1, nn);
    //contactStat[nn]++;
  }

  delete nContacts; delete nContactsR; delete distances; distances=NULL;

  //if (length_contact <=totalContact+2)
  //  cerr << "  ^^^ Got more contact than expected: " << length_contact << " : " << totalContact<< endl;
  nDataPoints=length_contact=totalContact; 
  if (length_contact==0)
     return -1;

  predicted_label = new int[nDataPoints];
  pairwise_features = new ScoreMatrix(1, nDataPoints, "pairwise_features");
  for ( int t = 0; t<length_contact; t++){
    int offset1 = centerMap[t], offset2 = rIdxMap[t];
    (*pairwise_features)(0,t) = (offset2 - offset1 + 0.000001) / length_seq;
  }
  return 0;
}

// Features: left-end vector([window_size/2]), features read from input 
// file (dim_one_pos*window_size), right-end vector ([window_size/2]), 
// sequence distance; sequence distance/length_seq.
Score ContactNumberSequence::getFeatures(int pos, int i)
{
  int pivot=pos*2;
  int offset1 = _features[pivot];
  int offset2 = _features[pivot+1];
  int iif=i-window_size/2*2;
  int idx = iif / m_pModel->dim_one_pos;
  Score ret=1; // when i==dim_features-3
  if (i<window_size/2*2){
    offset1+=window_size/2;
    offset2+=window_size/2;
    //if (proc_id==0) cerr << "offset1="<<offset1<<", offset2="<<offset2<<endl;
    if ((offset1<window_size/2 and i==offset1) or (offset2<window_size/2 and i==offset2))
      return 1;
    else if ((length_seq-offset1==window_size/2*2-i and offset1>=length_seq-window_size/2)
            or (length_seq-offset2==window_size/2*2-i and offset2>=length_seq-window_size/2))
    {
      return 1;
    } else 
      return 0;
  } else if (idx<window_size){
    if (offset1+idx<0 or offset1+idx>=length_seq) return 0; // training Secondary Structure: allow termial residues
    ret= obs_feature[offset1+idx][iif%m_pModel->dim_one_pos];
  } else if (i<dim_features-3){
    idx-=m_pModel->window_size;
    ret= obs_feature[offset2+idx][iif%m_pModel->dim_one_pos];
  } else if (i==dim_features-2)
    ret = Rg/10;
  else if (i==dim_features-1)
    ret = abs(offset2-offset1+0.000001)/length_seq;

  return ret;
}

Score ContactNumberSequence::MAP(){ // Posterior Decoding (Marginal Probability Decoder)
  memset(pred, 0, sizeof(int)*20);
  Score bestScore, obj = 0;
  Score *vi2;
  int i, t, bestState=0; 

  for(t=nDataPoints-1; t>=0;t--)
    obj+= (*arrVi)(t, obs_label[t]) - Partition[t];

  for(t=nDataPoints-1, vi2=&((*arrVi)(t, 0)); t>=0; t--, vi2-=arrVi->cols){
    bestState=0; 
    bestScore=vi2[0];
    for(i=1; i< num_states; i++){
      if (vi2[i]>bestScore){
        bestScore=vi2[i];
        bestState=i;
      }
    }
    predicted_label[t]=bestState;
    pred[predicted_label[t]]++;
  }
  cerr << " ** ** Label distribution: "; 
  for (int i=0; i<num_labels-1; i++) 
    cerr << pred[i] << ", "; 
  cerr << pred[num_labels] << endl;
  return obj;
}

void ContactNumberSequence::Predict()
{
  int i, j, t, startRes=m_pModel->startRes;
  int endRes=m_pModel->endRes;
  bool bGap;
  if (startRes<0)
    startRes=1;
  if (endRes<1)
    endRes=length_seq;

  contact_DIST=3;
  cerr << "  --------- " << p_name << "; len_contact=" << length_contact << endl;
  ComputeGates();
  ComputeVi();
  CalcPartition();
  MAP();

  for ( i=0; i<num_labels; i++)
    for ( j=0; j<num_labels; j++) {
      confusion[i][j]=0; confusionEnergy[i][j]=0;
    }

  // comparison
  int* num_contact8 = new int[length_seq];
  memset(num_contact8, 0, sizeof(int) * length_seq);
  int head = max(0,startRes-1), tail = min(length_seq, endRes);
  fprintf(stderr, "\n\tSeq :: (leftRes, cn)-(rightRes, cn) head=%d; tail=%d\n", head, tail);
  for (t=0; t<length_seq; t++){
    bGap=false;
    for (int g=0; g<(int)gaps.size(); g++)
      if ( abs(gaps[g]-t) <= window_size/2) { 
        bGap=true; break; 
      }
    int contactH=_MAX_CONTACT-1, contactT=_MAX_CONTACT-1;
    int idxH=0, idxT=length_seq-1;
    // cerr << " ===== " << t << " :: ";
    for (i=t-contact_DIST; i>=0; i--){
      int idx=IdxMap[t*length_seq+i];
      // cerr << idxH << ", " << predicted_label[idxH] << "," << contactH << "  ";
      if (idx>=0 and predicted_label[idx]<=contactH){
        contactH=predicted_label[idx];
        idxH=idx;
        //if (contactH<num_labels-1)
        //  cerr << "("<<i<<","<<num_labels-contactH-1 << ")->";
      }
    }
    //cerr << "\tor\t";
    for (i=t+contact_DIST; i<tail; i++){
      int idx=IdxMap[t*length_seq+i];
      if (idx>=0 and predicted_label[idx]<=contactT){
        contactT=predicted_label[idx];
        idxT=idx;
        // if (contactT<num_labels-1) 
        //   cerr << "("<<i<<","<<num_labels-contactT-1 << ")=>";
      }
    }
    // cerr << "\tpredictCN " << t << ":\t" << head << ":" << idxH << ":" 
    //      << contactH << "; " << tail << ":" <<idxT << ":" << contactT << endl;
    num_contact8[t]=num_labels*2-2-(contactH+contactT);
    //cerr << "\t::" << num_contact8[t] << " vs " << obs_cnum[t] << endl;
    for (int k=0; k<_MAX_CONTACT; k++) {
      (*prob_residue2state)(t,k) = 0;
      for (int k1=0, k2=k; k1<=k; k1++, k2--)
        (*prob_residue2state)(t,k)+=(*arrVi)(idxH, k1)*(*arrVi)(idxT, k2);
      (*log_prob)(t, k) = log((*prob_residue2state)(t,k));
    }
  }
  for (t=2; t<length_seq-2; t++){
    bGap=false;
    for (int g=0; g<(int)gaps.size(); g++)
      if (abs(gaps[g]-t)<=window_size/2){ bGap=true; break; }
      //if (gaps[g]==t){ bGap=true; break; }
    //if(bGap) continue;
    if (num_contact8[t] == 0)
      num_contact8[t] = (num_contact8[t-2] + num_contact8[t-1] + num_contact8[t+1] 
                                           + num_contact8[t+2]) / 4;
    else if (num_contact8[t] > 9)
      num_contact8[t] = 9;
    pred_cnum[t] = min(_MAX_CONTACT-1, num_contact8[t]);
  }
  FILE* fp=stdout;
  fprintf( fp, "\n>Prediction of Contact Number %s (length of sequence=%d:%d-%d)\n", 
           p_name.c_str(), length_seq, startRes, endRes);
  fprintf( fp, ">residueNo\t<8A\tTureNo.\n");
  
  double sumx=0, sumx2=0, sumy=0, sumxy=0, sumy2=0, sumR=0;
  int x, y, l=0;
  int len=0; //= length_seq-gaps.size()-contact_DIST*2;
  double square_err=0, ab_err=0;
  for (t=0; t<length_seq; t++){
    bGap=false;
    for (int g=0; g<(int)gaps.size(); g++)
      if (abs(gaps[g]-t)<=window_size/2){ bGap=true; break; }
      //if (gaps[g]==t){ bGap=true; break; }
    //if(bGap) continue;
    x=pred_cnum[t], y=obs_cnum[t];
    if (obs_exist_flag)
      fprintf(fp, ">%9d%9d%9d\n", t+1, x, y);
    else
      fprintf(fp, ">%9d%9d\n", t+1, x);

    square_err += (x-y)*(x-y);
    ab_err += abs(x-y);
    if (y>0){
      sumR += abs(x-y)/y;
      l++;
    }
    sumx+=x, sumx2+=x*x;
    sumy+=y; sumy2+=y*y;
    sumxy+=x*y;
    len++;
  }
  //double corr= (len*sumxy - sumx*sumy) / sqrt((len*sumx2-sumx*sumx)*(len*sumy2-sumy*sumy));
  //double DevA=sqrt((sumy2-2*sumxy+sumx2)/(sumy2-sumy*sumy/len));
  double DevR=sumR/l;
  l=len;
  double corr= (l*sumxy - sumx*sumy) / sqrt((l*sumx2-sumx*sumx)*(l*sumy2-sumy*sumy));
  double DevA=sqrt((sumy2-2*sumxy+sumx2)/(sumy2-sumy*sumy/l));
  if (obs_exist_flag)
    cout << p_name << "(" << len << ") sqrt(square_err)= " << sqrt(square_err/len)
         << "; mean(err)= " << ab_err/len << "; corr= " << corr << "; DevA= " << DevA 
         << "; DevR= " << DevR << endl;

  m_sumx=sumx, m_sumx2=sumx2, m_sumy=sumy, m_sumxy=sumxy, m_sumy2=sumy2, m_sumR=sumR, m_len=len;
  delete num_contact8;
}

double ContactNumberSequence::CalculateEnergy(int res_no, int cnum)
{
  int obs_label=min(_MAX_CONTACT-1, cnum);
  if (res_no>=length_seq or res_no<0)
    return 0;
  int pred_label=pred_cnum[res_no];
  confusion[pred_label][obs_label]++;
  //if (pred_label<_MAX_CONTACT) { pairs_rmsd++; rmsd_total+=abs(obs_label-pred_label); }

  double log_rf=logReferenceState[cnum];
  double energy= (*log_prob)(res_no, obs_label) - log_rf;
  double forbidden_val=10.0;
  if (energy>forbidden_val)
    energy=forbidden_val;
  
  return energy;
}

PNN1Inference::PNN1Inference()
{
}

PNN1Inference::~PNN1Inference()
{
  for (vector<SEQUENCE*>::iterator it=testData.begin(); it<testData.end(); it++)
    delete *it;
}

void PNN1Inference::Predict(string input, int Format)
{
  if (Format==0)
    LoadData(input);
  else
    LoadDataSVMFormat(input);

  totalPos = totalCorrect = 0;
  for(int i=0;i<num_labels;i++) 
    for(int j=0;j<num_labels;j++) 
      confusion[i][j]= confusion_sum[i][j]=0;

  for ( vector<SEQUENCE*>::iterator it = testData.begin(); it < testData.end(); it++){
    (*it)->Predict();
  }

  //testData->ComputeTestAccuracy();

  cout << " test ACC(MAP): " << (double) totalCorrect/totalPos 
       << "   " << totalCorrect << "/" << totalPos << endl;

  cout << endl << "Confusion Matrix: " <<endl;
  memset(sum, 0, sizeof(int)*30);
  for(int i=0;i<num_labels;i++){
    double tot=0;
    for(int j=0;j<num_labels;j++){
      printf("%7d ",confusion_sum[i][j]);
      tot+=confusion_sum[i][j];
      sum[j]+=confusion_sum[i][j];
    }
    printf("%7.2f\n", 100.0*confusion_sum[i][i]/tot);
  }
  for (int i=0; i<num_labels; i++)
    printf("%7.2f ", 100.0*confusion_sum[i][i]/(0.0001+sum[i]));
  cout << endl;
}

void PNN1Inference::SetParameters(int model, int n_local, int n_pairwise)
{
  this->model = (MODEL) model;
  if (this->model==singleton)
    dim_features = window_size/2*2 + window_size*dim_one_pos + 1; // Secondary Structure
  else if (this->model==pairwise)
    dim_features = window_size/2*2 + (2*window_size*dim_one_pos) + n_pairwise + 1; // constant 1 
  else
    // 2+1 means Rg, chain distance(cd)/length_seq, and constance 1 
    dim_features = window_size/2*2 + (2*window_size*dim_one_pos) + 2 + 1; 
  
  num_states=num_labels;
  iGateStart = 0;
  iFeatureStart = iGateStart+num_states*num_gates;
}

// return the total number of valid columns
int PNN1Inference::getColumns(string columns, int* column_flags)
{
  memset(column_flags, 0, sizeof(int)*MAX_COL);
  if (columns.length()>0) {
    //cerr << "|" << columns << "|" << endl;
    istringstream si(columns);
    char c=' ';
    int col1,col2;
    while (!si.eof()){
      si >> col1 >> c;
      //cerr << "**" << col1 << ":" << c << endl;
      if (c==','){
        if (col1>=0 and col1<MAX_COL) column_flags[col1]=1;
      } else if (c=='-'){
        si >> col2 >> c;
        if (col1>=0 and col1<MAX_COL){
          for (int i=col1; i<=min(MAX_COL-1, col2); i++)
            column_flags[i]=1;
        }
      } else {
        if (col1>=0 and col1<MAX_COL) column_flags[col1]=1;
        break;
      }
    }
    //cerr << "selected columns: ";
    //for (int i=0; i<MAX_COL; i++)
    //  if (column_flag[i]==1) cerr << i << ",";
    //cerr << endl;
  }
  int total_col=0;
  for (int i=0; i<MAX_COL; i++)
    total_col+=column_flags[i];
  return total_col;
}

int PNN1Inference::Initialize(int model, string model_dir, string columns,
                              string pairwise_columns, string init_file)
{
  model_file = model_dir+init_file;

  cerr <<"num_params = " << num_params << endl;
  cerr << "  init_file: " << init_file << endl;
  ifstream fin(init_file.c_str()); // no file header
  char buf[10240];
  string tmp;
  int data;
  for (int i=0; i<7; i++) {
    fin.getline(buf,10240);
    if (fin.gcount()==0)
      break;
    else {
      istringstream si(buf);
      si >> tmp >> data;
      if (tmp=="num_params:")
        num_params = data;
      else if (tmp=="num_gates:")
        num_gates = data;
      else if (tmp=="window_size:")
        window_size = data;
      else if (tmp=="dim_features:")
        dim_features = data;
      else if (tmp=="num_labels:")
        num_labels = data;
      else if (tmp=="dim_one_pos:")
        dim_one_pos = data;
      cerr << tmp << "::::"<<data<<endl;
    }
  }
  weights = new double[num_params];
  cerr << strlen(buf) << "  Last 10 weights are : " << endl << "    ";
  for(int i=0;i<num_params;i++) {
    fin >> weights[i];
    if (i>num_params-10) cerr << weights[i] << ", ";
  }
  fin.close();
  cerr << endl;

  int n_local = 1 + getColumns(columns, column_flags);
  int n_pairwise = getColumns(pairwise_columns, pairwise_column_flags);
  SetParameters(model, n_local, n_pairwise);

  return 0;
}

int PNN1Inference::LoadData(string input)
{
  totalPos = totalCorrect = 0;
  for(int i=0;i<num_labels;i++) 
    for(int j=0;j<num_labels;j++) 
      confusion[i][j]= confusion_sum[i][j]=0;

  map<string, int> pmap; // list of proteins that we need to calculate
  char buf[101];
  ifstream listf("../Data/contactNumberCompare.lst");
  while (!listf.eof()){
    listf.getline(buf, 100);
    string pname = buf;
    pmap[pname.substr(0,4)]=1;
    //cerr << pname << ":" << pmap.find(pname)->second << endl;
  }

  vector<int> myseq_test;

  // Read the training data
  ifstream trn_in(input.c_str());
  trn_in >> num_data;
  cerr << input << ": gateStart=" << iGateStart << "; ifs=" << iFeatureStart 
       << " - Load Data: "<< dim_one_pos << ";" << num_data << "; model=" << model 
       << "; pairwise=" << dim_pairwise_features << endl;
  
  double sumx=0, sumx2=0, sumy=0, sumxy=0, sumy2=0, sumR=0;
  int len=0;

  bool title_line_read=false;
  double dist, tmp;
  int a1, a2, length_s, length_labels, real_dim_one_pos, gap;
  string stmp, sline;
  SEQUENCE *seq = NULL;
  getline(trn_in, sline);
  for (int i=0; i<num_data; i++) {
    if (!title_line_read)
      getline(trn_in, sline);
    else
      title_line_read = false;
    istringstream iss(sline);
    iss >> length_s >> length_labels>> stmp;
    // construct a new sequence
    switch (model)
    {
    case singleton:  // e.g. Secondary Structure
      seq = new SingletonSequence(length_s, this, stmp);
      break;
    case pairwise:
      seq = new PairwiseSequence(length_s, this, stmp);
      break;
    case DistancePotential:
      seq = new DistancePotentialSequence(length_s, this, stmp);
      break;
    case contact_number:
      seq = new ContactNumberSequence(length_s, this, stmp);
      break;
    }
    // cerr << i<<":"<<length_s<<" " << stmp <<endl;

    for(int j=0;j<length_s;j++){
      getline(trn_in, sline); //cerr << "Contact " << j << ":";
      istringstream iss(sline);
      int d=0;
      for (int k=0; !iss.eof(); k++) {
        iss >> tmp;
        if (column_flags[k]==1) {
          seq->obs_feature[j][d++] = tmp; //cerr << tmp <<" ";
        }
      } //cerr << endl;
      real_dim_one_pos=d+1;
      seq->obs_feature[j][d] = 1;
    }
    cerr << " Data of residue " << i <<endl;
    
    if (length_labels<=0)
      seq->obs_exist_flag=false;

    if (model==singleton){
      for(int j=0;j<length_labels;j++){
        getline(trn_in, stmp); //cerr << "Contact " << j << ":";
        if (stmp[0]=='#'){
          gap=atoi(stmp.substr(5).c_str());
          seq->gaps.push_back(gap);
        } else {
          istringstream iss(stmp);
          iss >> a1; // if (!proc_id) cerr << a1 << "(" << j << ")" << endl;
          seq->obs_label[j]=a1;
        }
      } //exit(0);
    } else { // pairwise
      getline(trn_in, stmp);
      for(int j=0;j<length_labels;j++){
        getline(trn_in, stmp); //cerr << "Contact " << j << ":";
        if (*(stmp.end()-1)=='#') {
          sline = stmp;
          title_line_read=true;
          break;
        }
        if (stmp[0]=='#'){
          gap=atoi(stmp.substr(5).c_str());
          seq->gaps.push_back(gap);
        } else if (stmp.size()>0){
          istringstream iss(stmp);
          iss >> dist >> a1 >> a2; 
          try {
            seq->distances->helper=dist;
            if (a1<length_s and a2<length_s){
              (*seq->distances)(a1,a2)=dist;
              (*seq->distances)(a2,a1)=dist;
            } else 
              cerr << stmp << ":: " << seq->p_name << " dist" 
                   << "(" << a1 << "," << a2 << ")=" << dist << endl;
            //cerr << (*seq->distances)(a1,a2) << "(" << a1 << "," << a2 << ") " << seq->p_name << endl;
            for (int k=0, d=0; !iss.eof(); k++) {
              iss >> tmp;
              if (pairwise_column_flags[k]==1) {
                (*seq->obs_pairwise_features)(a1, a2, d++) = tmp; //cerr << tmp <<" ";
              }
            } //cerr << endl;
          } catch (bad_exception e){
            cerr << seq->p_name << " dist" << "(" << a1 << "," << a2 << ")=" << dist << endl;
            exit(0);
          }
        }
      }
    }
    //string p_name = seq->p_name.substr(0, 4);
    //p_name[p_name.size()-1]=tolower(p_name[p_name.size()-1]);
    //if (pmap.count(p_name)==0) {
    //  cerr << " Ignore " << p_name << " " << i << "; length=" << seq->length_seq << endl;
    //  delete seq; 
    //  continue;
    //}
    cerr << " Insert " << i << "; length=" << seq->length_seq << endl;
    //testData.push_back(seq);
    if (seq->makeFeatures()!=-1){
      seq->Predict();
      sumx += seq->m_sumx, sumx2 += seq->m_sumx2, sumy += seq->m_sumy,
                           sumxy += seq->m_sumxy, sumy2 += seq->m_sumy2,
                           sumR += seq->m_sumR, len += seq->m_len;
    }
    delete seq;
  }
  trn_in.close();

  double DevR=sumR/len;
  int l=len;
  double corr= (l*sumxy - sumx*sumy) 
             / sqrt( (l*sumx2-sumx*sumx) * (l*sumy2-sumy*sumy) );
  double DevA=sqrt( ( sumy2 - 2*sumxy+sumx2) / (sumy2-sumy * sumy / l) );
  cout <<  " Total Corr/devA/DevR..(" << len << ") ; corr= " << corr 
       << "; DevA= " << DevA << "; DevR= " << DevR << endl;

/*
  if (model==contact_number){
    cerr << "\ttest Data" << endl;
    for (int i=0; i<num_labels; i++)
      cerr << i << ": " << contactStat_sum[i] << endl;
    //exit(0);
  }
*/
}

void PNN1Inference::LoadDataSVMFormat(string input)
{
  if (model!=singleton)
    return; // svm format only works for singletons

  totalPos = totalCorrect = 0;
  for(int i=0;i<num_labels;i++) 
    for(int j=0;j<num_labels;j++) 
      confusion[i][j]= confusion_sum[i][j]=0;

  // Read the training data
  int total_data = 0;
  string sline, token;
  vector<string> buffer;
  ifstream trn_in(input.c_str());
  while (!trn_in.eof()) {
    getline(trn_in, sline);
    if (sline.size()==0)
      continue;
    buffer.push_back(sline);
    total_data ++;
  }
  num_data = 1;
  double tmp;
  int length_s=total_data/num_data, real_dim_one_pos;
  cerr << input << ": gateStart=" << iGateStart << "; ifs=" << iFeatureStart 
       << " - Load Data: " << dim_one_pos << ";" << num_data << "; model="
       << model << "; " << total_data << endl;
  SEQUENCE *seq;
  for(int i=0;i<num_data;i++){
    // construct a new sequence
    seq = new SingletonSequence(length_s, this, "");
    for(int j=0;j<length_s;j++){
      sline = buffer[i*length_s+j];
      int lastspace=0, lastcolon=0, idx=0, k=0, labeled=0;
      for ( ; k<sline.size(); k++) {
        if (sline[k]==' '){
          if (!labeled) {
            string token=sline.substr(lastspace, k-lastspace);
            seq->obs_label[j]=atoi(token.c_str())-1;
            if (j==length_s-1) cerr<<endl<<j<<" "<<seq->obs_label[j]<<" ";
            labeled=1;
          } else {
            if (lastcolon==0){
              idx++;
              token=sline.substr(lastspace, k-lastspace);
            } else 
              token=sline.substr(lastcolon, k-lastcolon);
            seq->obs_feature[j][idx-1] = atof(token.c_str()); 
            //if (j==0) cerr<<idx<<":"<< seq->obs_feature[j][idx-1] <<" ";
          }
          lastspace = k+1;
        } else if (sline[k]==':') {
          string token=sline.substr(lastspace, k-lastspace);
          idx = atoi(token.c_str());
          lastcolon = k+1;
        }
      }
      if (lastcolon==0){
        idx++;
        token=sline.substr(lastspace, k-lastspace);
      } else 
        token=sline.substr(lastcolon, k-lastcolon);
      seq->obs_feature[j][idx-1] = atof(token.c_str()); 
      //if (j==0) cerr<<idx<<":"<<seq->obs_feature[j][idx-1] <<" ";
      real_dim_one_pos=idx+1;
      seq->obs_feature[j][idx] = 1;
    }
    cerr << endl << " Predicting " << i << " length=" << seq->length_seq << endl;
    if (seq->makeFeatures()!=-1)
      seq->Predict();
    delete seq;
  }
  dim_one_pos = real_dim_one_pos;
  trn_in.close();
}

void Usage()
{
  cerr << "PNN1Inference [-i input_file] [-d model_dir] [-m model] "
       << "[-l columncs] [-p pairwise_features] [-f model_file] [-F format]\n";
  cerr << "\t-m model: singleton, pairwise, DistancePotential or contact_number" << endl;
  cerr << "\t-l columns: select those columns to calculate in the singleton "
       << "features section. Note: the columns start from 0. For example, "
       << "-l0,3-6 means columns 0, 3, 4, 5 and 6." << endl;
  cerr << "\t-p pairwise_feaures : select those columns to calculate in the "
       << "pairwise feature section. Similar to -l option. Default value "
       << "is -1, meaning no feature." << endl;
  cerr << "\t-F format: 0 denotes sequence format; 1 denotes svm format. Default is 0." << endl;
  cerr << "  Example: " << endl;
  cerr << "     ./PNN1Inference -i ss70 -d . -l 0-39 -m singleton -f model100" << endl;
  cerr << "  This command runs a 3 state secondary structure classfier. " << endl;
  cerr << "  As shown in model100 file, the window size is set to be 5, "
       << "using 20 gates, 40 features for each residue, i.e. column 0 to 39." << endl;
  cerr << "  ss70 file has 513 proteins sequences. We use it as the test set." << endl;
}

int main(int argc, char **argv){

  // the command line must be:
  if (argc < 3) {
    Usage();
    return 0;
  }

  extern char* optarg;
  char c=0;
  int model=1, Format=0;
  string model_dir="./", input_f="ss70", columns="", pairwise_columns="";
  string init_file="";
  while((c=getopt(argc,argv,":d:i:l:p:f:m:F:"))!=EOF){
    switch(c)
    {
    case 'i':
      input_f=optarg;
      break;
    case 'l':
      columns=optarg;
      break;
    case 'p':
      pairwise_columns=optarg;
      break;
    case 'd':
      model_dir = optarg;
      if (model_dir[model_dir.size()-1] != '/')
        model_dir += "/";
      //cerr << "model_dir=" << model_dir << endl;
      break;
    case 'f':
      init_file = optarg;
      break;
    case 'm':
      if (strcasecmp(optarg, "singleton")==0)
        model = 0;
      else if (strcasecmp(optarg, "pairwise")==0)
        model = 1;
      else if (strcasecmp(optarg, "DistancePotential")==0)
        model = 2;
      else if (strcasecmp(optarg, "contact_number")==0)
        model = 3;
      break;
    case 'F':
      Format = atoi(optarg);
      //cerr << " Format=" << Format << "=" << optarg << endl;
      break;
    default:
      Usage();
      exit(-1);
    }
  }

  PNN1Inference cnfModel;
  cnfModel.Initialize(model, model_dir, columns, pairwise_columns, init_file);
  //cerr << "Initialization Finished! " << endl;
  //cerr << "Format is " << Format << endl;
  //cnfModel.Predict(input_f, Format);
  if (Format==0)
    cnfModel.LoadData(input_f);
  else
    cnfModel.LoadDataSVMFormat(input_f);
  return 0;
}
